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ABSTRACT 


IMPLEMENTATION OF UNIFORM PERTURBATION METHOD FOR POTENTIAL 
FLOW PAST AXISYMMETRIC AND TWO-DIMENSIONAL BODIES 

Tin Chee Wong 

Old Dominion University, 1984 
Director: Dr. Surendra N. Tiwari 
Dr. Chen-Huei Liu 

The aerodynamic characteristics of potential flow past an 
axi symmetric slender body and a thin airfoil are calculated by using a 
uniform perturbation analysis method. The method is based on the 
superposition of potentials of point singularities distributed inside 
the body. The strength distribution satisfies a linear integral 
equation by enforcing the flow tangency condition on the surface of the 
body. The complete uniform asymptotic expansion of its solution is 
obtained with respect to the slenderness ratio by modifying and adapting 
an existing technique. Results calculated by the perturbation analysis 
method are compared with the existing surface-singularity -panel method 
and some available analytical solutions for a number of cases under 
identical conditions. From these comparisons, it is found that the 
perturbation analysis method can provide quite accurate results for 
bodies with small slenderness ratio. The present method is much simpler 
and requires less memory and computation time than an existing surface- 
singularity-panel method of comparable accuracy. 
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Chapter 1 
INTRODUCTION 

In recent years, applications of powerful and sophisticated 
computers are being made extensively to analyze and solve important 
aerodynamic problems. Computer codes are now being developed to 
simulate the flowfield over entire aircraft configurations by 
considering the boundary -layer and full -potential formulations. For 
realistic physical applications. It would be desirable to work with the 
full Navier-Stokes equations. However, this requires large computer 
resources (memory and execution time) which makes the computational 
costs prohibitive for many applications. For this reason, the former 
method is highly desirable for solving many important aerodynamic 
problems. Thus, potential solutions become essential in analyzing many 
flow problems. Although it would be desirable to obtain potential 
solutions for the complete aerodynamic configurations, it is usually 
convenient to analyze different basic components of the system 
separately. In analyzing the flowfield over an airplane, it is 
customary to divide the entire problem into the three separate problems 
of fuselage, wing, and tail. Thus, separate programs are developed to 
analyze the flow over the fuselage, wing and tail. The analyses of 
axl symmetric bodies and two-dimensional airfoils are considered as the 
first step for fuselage and wing designs respectively. The two- 
dimensional airfoil results are modified usually to account for the 
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effect of three-dimensional wing flow. 

The incompressible potential flow over axi symmetric bodies and 
wings is essentially governed by the Laplace equation. Although the 
Laplace equation is one of the simplest and well known of all partial 
differential equations, the number of known closed-form solutions is 
quite small. This is due to the non -homogeneous boundary conditions for 
the partial differential equation and the requirement for the body 
surface to be a coordinate surface of the one of the orthogonal 
coordinate system for which the Laplace equation can be solved by 
separation of variables. As a consequence, approximate solutions become 
more useful in practical applications. One typical example is the 
slender-body theory where the body boundary condition is simplified for 
small slenderness ratio and angle of attack. In many Important fields 
successful approximations are found and some of these are discussed in 
[1]*. The present study is an application of the slender body theory. 

Many previous Investigators have attempted to solve the potential 
flow problem around a slender body. The idea of using a distribution of 
singularities Inside the body was first treated by von Karman [2]. For 
the axi symmetric motion, he used a continuous distribution of sources 
along the axis of the given body and gave a method for computing the 
distribution strength. For the lateral flow past a given body he used a 
continuous distribution of doublets along the axis with the doublet axes 
normal to the body axis. It Is difficult to obtain a closed-form 
solution for the exact problem; however, an approximate closed-form 
solution may be obtained, if the shape of the body satisfies certain 


♦The numbers in brackets indicate references 
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conditions. The body should be slender and should not have any 
discontinuities in its surface slope. Zedan and Dalton [3,4] and 
Kuhlman and Shu [5,6] extended the method by dividing the singularity 
distribution into elements and employing piecewise linear or polynomial 
functions to represent the variations of the Intensity of the source 
distribution over each line element. The effects of the order of the 
distribution, the number of elements, the normalization of the body 
coordinates, the slenderness ratio and the geometry of the profile on 
the performance of the method were studied in detail. It was concluded 
that each of these factors should be chosen appropriately in order to 
have accurate results. 

Another approach to obtain approximate solutions to the slender- 
body problem is to use the perturbation method. In 1967 Handel sman and 
Keller [7] introduced a method to give a uniformly-valid approximate 
solution of axially symmetric potential flow containing a small 
parameter, known as slenderness ratio e. The source strength distri- 
bution satisfied a linear integral equation. A uniform asymptotic 

solution to the Integral equation was obtained. The left-hand side of 

2 

the integral equation was expanded In the power series of e and the 

2 2 
right-hand side in the power series of e and a power series in e 

2 2 

multiplied by e (log e ) without taking account of the dependence of the 
source strength on e. It was assumed also that the source strength 
vanishes in an Interval near each end of the body and that the profile 
curve of the body is analytic. It was possible to give a complete 
uniformly-valid solution to any order of approximation in this way. The 
method used to determine the axially symmetric potential flow about a 
slender body of revolution is presented in [7]. 
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The axial singularity method can be applied to a two-dimensional 
airfoil in a similar manner to solve the potential flow problem. 
Similarly, one can have an alternative approach that will simplify the 
mathematical conditions of the problem and enable the construction of an 
approximate closed-form solution. Such an approach is called the 
classical thin airfoil theory. The thin airfoil theory breaks down near 
stagnation points and near the edges of the airfoil. This is because 
the assumption of small perturbation is not valid in the neighborhood of 
such points; the perturbation velocity is of the same order of magnitude 
as the undisturbed velocity. For a survey of these difficulties and the 
methods that are developed to relieve them one should make reference to 
Van Dyke [1]. The first approach is the method of matched asymptotic 
expansions and the second is the method of strained coordinates [8]. 

The method of strained coordinates has been modified by Hoogstraten 

[9] . However, these methods are so complicated that only the first two 
terms of a uniform expansion are generally obtained. Geer and Keller 

[10] presented the method to obtain a uniform asymptotic solution for 
potential flow around a thin symmetric airfoil by adopting the method of 
Handel sman and Keller [7]. This Idea was extended and modified by Geer 

[11] to obtain the complete uniform asymptotic expansion for the non- 
symmetric airfoil. This perturbation method Is presented in Chap. 2. 
Since the method to obtain the complete uniform asymptotic expansion for 
potential flow past axi symmetric and two-dimensional bodies is known, 
the only interest here is the assessment of the perturbation method. 

Such questions as: how large can the slenderness ratio be or what are 

the advantages and disadvantages of this method over other existing 
methods are addressed. In the light of these considerations, some 
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programs have been developed to obtain the pressure distributions for 
slender axi symmetric bodies and thin airfoils. In order to make 
comparisons, some exact solutions and solutions obtained by using an 
existing surface-singularity-panel method have been formulated and 
presented in the latter part of Chap. 2. Much effort has been devoted 
to obtain the expressions of the asymptotic expansions for the potential 
function, with the emphasis on some axi symmetric slender bodies and thin 
airfoils, and these are presented in Chap. 3. The pressure distrbutions 
on axi symmetric bodies and airfoils obtained by the perturbation 
analysis method are compared with the results generated by the panel 
method as well as exact solutions whenever available. These results are 
presented in Chap. 4. The perturbation analysis method can be applied 
to solve the electrostatic and magnetostatic potentials around a slender 
conducting body [10-13]. Recently, Homentcovschi ever* modified this 
method to solve two-dimensional elasticity problems [14]. 
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Chapter 2 

THEORETICAL FORMULATION 

Basic potential flow theory is reviewed in this chapter. The 
problem under consideration in this study is that of the steady 
irrotational flow of incompressible, inviscld fluid, for which the 
Navier-Stokes eqautions can be reduced to the classical Laplace 
equation. Consider a body Immersed in an ideal fluid and a space-fixed 
reference frame. Then the surface of the body is described by 

F(x 1 ,x 2 ,x 3 ) = 0 (2.1) 

The variables x ^ ,x 2 and x 3 are the Cartesian coordinate system used in 
the analysis. Because the fluid is incompressible, the continuity 
equation implies that the divergence of velocity must be zero 

V.V * 0 (2.2) 

Helmholtz's theorem states that the vortlcity of an initially 
irrotational, invlscid fluid is zero. Mathematically, one has 

VxV * 0 (2.3) 


and it follows that 



7 


V = v+ (2.4) 

Here <f> Is the velocity potential. The substitution of Eq. (2.4) into 
Eq. (2.2), gives the governing Laplace differential equation for 
incompressible potential flow, i.e.. 


V 2 * = 0 (2.5) 

In order to solve for the potential 4 , Eq. (2.5) should satisfy the 
boundary conditions 


V$. VF ■ 0 on F(xj,X 2 »X 3 ) = 0 (2.6) 


Also, the components of the velocity should vanish in a certain way with 
distance from the body as that distance tends to infinity [15]. The 
Bernoulli equation relating pressure and velocity is given by 



P - P 


T 

7 


pU 


= 1 


U 2 

- <TT> 


(2.7) 


Once the velocity on the surface of the body, U s , is determined, the 
pressure distribution can be computed through Eq. (2.7). 


2.1 Perturbation Analysis Technique 
A large number of perturbation techniques to obtain approximate 
solutions to physical potential flow problems is available in the 
literature [1,7-12]. In this chapter, the special perturbation 
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technique developed by Handel sman and Keller [7,10-12] is reviewed in 
detai 1 . 


2.1.1 Axially Distributed Singularity Method for Axisymmetric Body 

Consider the steady state of incompressible, inviscid and 
irrotational fluid past an axially symmetric rigid body which is at zero 
angle of attack to the uniform stream, see Fig. 2.1. It is convenient 
to analyze the problem by introducing cylindrical coordinates (r,e,x) 
with origin at the body's nose and the x-axis along its line of 
symmetry. The surface of the body is described by an equation r 
= e R(x), where e is the slenderness ratio, i.e., the ratio of the 
maximum radius of the body to its length. If one uses the length of the 
body as a unit of length, then the body intercepts the x-axis at x=0 and 
x»l. Then, the maximum value of R(x) in absolute value is one. In this 
chapter, the slender body of revolution is considered, and e is 
considered to be a small value. 

Since it is assumed that the flow is irrotational , the velocity 
potential surrounding an axially symmetric body is given as 

♦ - «j>° + 4> b (2.8) 

where <j>° is the potential of the Incident flow while <(> b is the 
perturbation potential. The perturbation potential <|> b is due to the 
presence of the body, and this is to be determined. This is represented 
by a superposition of point sources distributed along the x-axis inside 
the body. The boundary conditions for this problem are the no- 
penetration (or tangent -flow) condition: 
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V<|> . n » 0 on r = eR(x) (2.9) 

and the Infinity condition: 

^ 0 as r infinity (2.10) 

where n is the unit outward normal vector to the surface of the body. 

Due to the axial symmetry of the flow, both <j> b and 4>° are independent 

of 9 and the equation of the profile curve depends on e, so that ♦ b 1sa 

function of x, r 2 and e. In this study, our aim is to obtain an 

b 2 

asymptotic expansion of <j» with respect to e around e * 0. We define 
S(x)=R 2 (x) and then the cross-sectional area of the body at x is 
ire 2 S(x). Assume that S(x) is analytic on 0 < x < 1, with S(0) = 0 = 

S(l) and it can be expanded in power series about endpoints as follows: 

S(x) - I c n x" (2.11) 

n=l n 

S(x) - I d (l-x) n (2.12) 

M wl ® * 


where c n = S^ n ^(0) /n i and d n * (-l) n S^ n ^(l) /n! . It is assumed that 
Cj and dj are non-zero, i.e., the radius of curvature at each end of the 
body is not zero. The perturbation potential ^ is used to represent 

the superposition of point sources distributed along a segment of the x- 
axis inside the body with unknown source strength density f(x,e) per 
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unit length [7]. Thus, 


<j>(x,r 2 ;e) = 4»°(x,r 2 ) - / 


8(e) f(S.e) d£ 

<x(e) C(x-C) 2 + r 2 ] 


172 


(2.13) 


where «(e) and 8(e) are constants which determine the extent of the 
source distribution and must be found in addition to f(x,e), and they 
satisfy the inequalities 0 < <* < 8 < 1. 

The potential function $ can be related to stream function i|> 
by i» x = -r<|> r and = r^. Thus, one can rewrite Eq. (2.13) as 



8(e) tx-e) ds 

ct(e) [(x-e) 2 + r 2 ] i/2 


(2.14) 


The boundary condition, Eq. (2.9), implies that there is no flow of 
fluid through the surface of the body, i.e.. 


8(e) 

/ f(S,e)d5 - 0 (2.15) 

o(e) 


Since the body is a continuation of the axial streamline for the total 
flow, one can obtain 


i»[x,e 2 S(x) ;e ] » 0 


(2.16) 


Using Eqs. (2.15) and (2.16) in Eq. (2.14), it follows that 


o, 2„, . _ 

+ [x,e S(x)] 


i f Bte) (x-e) f(E,s) de 
^ «(c) C(x-e) 2 + e 2 S(x)] 1/2 


(2.17) 


Equation (2.17) is a linear integral equation from which f(t;,e), o(e) 
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and 0(e) need to be determined. The uniform asymptotic expansion of 

the solution of this equation is obtained by adopting the method of 

Handel sman and Keller [7]. Since the left side of Eq. (2.17) is 

analytic in x for 0 < x < 1, f(x,e) must be analytic in its domain of 

definition a < x < 0, and the coefficients in the expansion of f(x,e) 

2 

with respect to e are assumed to be continuous in the Interval 

2 

0 < x < 1. Both o and 0 are found to be power series in e of the form 

a(e) = y a e 2n (2.18a) 

n=l n 

0(e) =1-1 0 e 2n (2.18b) 

n=l n 

In order to obtain an asymptotic expansion of the source strength 

2 

density f(x,e) of Eq. (2.17) with respect to e , one first expands each 

2 

side of Eq. (2.17) with respect to e without taking account of the 

2 

dependency of f on e . The left side can be expanded in a power series 

2 o p 

in e because ♦ is analytic In r . The right side can be expanded 

2 2 

asymptotically in powers of e and powers of e multiplied by 
2 

log(e ) with coefficients which are linear expressions in f. When both 
sides are expanded, Eq. (2.17) becomes 

00 j 04 X 0(e) 

4ir l )Ax)S (xje 23 = / f d£ - / f d5 
j=l J a(e) x 

*1 e y (L. t log e 2 G f )f 
j*l J J 


(2.19a) 
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where ^.(x) 3 jr (~^) J *°(x»r 2 ) | 2 


3r“ r = 0 


(2.19b) 


In Eq. (2.19a) Lj and Gj are linear operators defined by Eqs. (4.13)- 
(4.15) in [7], with e replaced by e 2 . Upon differentiating Eq. (2.19a) 
with respect to x, and by noting that 


d 

Hx 


[/ x 


0 


f dc - / f d$ = 2 f(x,e) 




( 2 . 20 ) 


there is obtained 


4. l 37 C+.CxJS^Cx)] - 2f(x,e) 

4,1 ox J 


+ A J ,f (2,n) 


To solve Eq. (2.21), the asymptotic solution for f is supposed to be of 
the form 


f(x,e) - I "f o 2 " (log t 2 ) m f_(x) 
n=l m=0 nm 


( 2 . 22 ) 


where f nm (x) are functions of x and must be determined by substituting 
Eq. (2.22) into Eq. (2.21), one obtains 


4 * l W - 2 i "i* * 2n 0og e 2 r f nm (x) 

n=i n=i m-u 


M 'i i E 2 ( j-4m)( log e 2)Hid (|_ , log e 2 e )f (x) 

n-1 m=0 j=l a * J J 


(2.23) 
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For each n and m , the coefficients of e 2n (log e 2 ) m on both sides of Eq. 
(2.23) can be obtained as follows: 


f 10 = 2 * 


(2.24a) 



n-j,0 


) 


n > 2 


(2.24b) 


f 



m > 1 


(2.24c) 


f nm = ' \ n>2 > m>1 (2 ‘ 2M) 

V. ■** J 


The f nm (x) can be determined recursively by starting with n=0, m=0. 
Thus, once the Lj and Gj are evaluated, Eqs. (2.22)-(2.24d) will yield 
the desired asymptotic solution for f(x,e) In terms of coefficients i|» n< 
Inserting Eq. (2.22) into Eq. (2.13), the asymptotic expansion of $ is 
obtained as 


♦(x,r 2 ;e) = <j>°(x,r 2 ) 


1 * n-1 n n 8 ( E ) 

- h l I £ Hog e 2 )" 1 / 

T n=l m=0 




.(.) [<*- 0 * 77 ?* 


(2.25) 


The method to obtain the asymptotic expansion with respect to e 
about e = 0 of the integral operator in Eq. (2.17) is available in 
C7]. 

Let the integral operator applied to a function F(x) which is 


Independent of e be 
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I(x,e) 


, 0 ( e) F 

«(«) C(x-5) 


.j -£)lx- 0 , d fi. 
2 + A(x)] 


17? 


(2.26) 


where a(e) and 0(e) are given In Eqs. (2.18a) and (2.18b). By adding 
and subtracting integrals to Eq. (2.26), one finds 


x B x 

I(x,e) - / F( 5 )dC - / F(|)dS + J F(e) 

a X a 




i&Si 
2 . 2 , 


iW 


-i]de 


b 

+ / F(0 

x 


[ 


lxj-.ll 


C(x-er+ e S(x)] 


C(x-5) + £ S(x)] 

+ 1 Me (2.27) 


']■ 


On the right-hand side of Eq. (2.27), one sets x-e»v and j--x=v to the 
third and fourth integral respectively. Then, Eq. (2.27) becomes 


I(x,e) * J* F(c)d5 - / 6(e) F(e)d? + W(x,e) + V(x,e) (2.28) 

a(e) x 

Here, W(x,e) = J F(x-v) [v(v 2 + e 2 S) -1 ^ 2 - 1] dv (2.28a) 

0 

B — X 

V(x,e) - - J F(x+v) [v(v 2 + e 2 S)“ 1/2 - 1] dv (2.28b) 

0 


To find the asymptotic expansion of W and V, the binomial expansion is 
considered 


y(y 2 + e 2 S) _1/2 - 1 - J a, (2.28c) 

j-1 v c 3 

where a^ * (-1) J (’2") (■§■) + j ” l)/jl • This expansion Is not 

valid throughout the domains of integration because these domains extend 
to v=0. Therefore, special treatment should be applied. To F(x+v) in 
the integrand of Eq. (2.28b), one adds and subtracts the two leading 
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2 2 - 1/2 

terms in its Taylor series at v=0. To [v(v + e S) ' - 1], one again 
adds and subtracts the leading terms of its binomial expansion. Thus, 
Eq. (2.28b) can be written as 


-V(x,e) = F(x)/ 0 * X [v(v 2 + £ 2 S )' 1/2 -l]dv + F'(x)/ M v[v(v 2 +e 2 S)' 1/2 -l]dv 
0 0 


+ I 


r B_X Cv(v - \ a^UCFU+v)- I F 

0 j=0 J j=0 


(j) 


+ a,e Z S / [F(x+v) - l F W v ] v~ z dv 

1 0 j*0 


(2.29) 


An analysis of the order of accuracy is described in Appendix A of 

2 2 

[7]. It shows that the first two integrals are 0(e) and 0(e log e ), 

respectively. The third term is 0(e ), but the fourth integral 

is 0(e ), so that it is asymptotically negligible compared to the three 

terms. To obtain the asymptotic expansion of the fourth integral in Eq. 

(2.29), one can apply the same procedure. To the first factor of the 

integrand, one adds and subtracts the second term in the binomial 

expansion given by Eq. (2.28c). To the second factor one can add and 

subtract the next two terms of the Taylor expansion of F(x+v) about 

v=0. In this manner, a sum of integrals is obtained; some of which can 

be evaluated explicitly, others can be expanded as power series 
2 

in e , and the remaining term is asymptotically small compared to the 
other terms. The same procedure can be applied to the remainder, and 
this can be 1 done repeatedly. Thus, V(x,e) can be expressed as a series 
of integrals which are successively smaller in e 2 a s e tends to zero. 
Furthermore, W(x,e) can be treated in an exactly analogous way, so one 
obtains 
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W(x,e) + V(x,e) = J a S n (x)e 2n [H (x,e) - H (x,e)] 

n=l n n n 

“ c(2n)/ \ ~ 

+ TSir'V" 1 - p n (x - e)] 

« c(2n+l), x 

- n I 0 CK n^ x * e) + 


(2.30) 


where, H n 


H n (x,e) = 


P n (x,£) = 
K n (x,e) = 


,P n f K n are defined by 


x “Ct « 2n-l r 

J v“ 2n [F(x-v) - l F 


0 

x-a 


j*l 


,2n 2j 2-% -1/2 


J v‘" [v(v‘+ £ ‘S) 
0 


J v 2n+1 [v(v 2 +e 2 S) -1 ^ 2 
o 



( a o =1 ) 


(2.30a) 

(2.30b) 

(2.30c) 


The functions H n> P n , and K n are defined by Eqs. (2.30a), (2.30b) and 

(2.30c) respectively, with (x-a) replaced by (8-x). To determine the 

asymptotic expansion of I(x,e), one must expand asymptotically the 

integrals Eqs. (2.30a) and (2.30b ) . H n (x,e) and H n (x,e) can be expanded 

o 

directly as Taylor series in e about e * 0. The integrals Eqs. (2.30b) 
and (2.30c) have been analyzed in detail in Appendix B of [7]. Appendix 
A in [7] shows that P n = 0(e 2n+1 ), K n = 0(e 2n+2 log e 2 ) and the same 

W «W fSJ 

holds true for P„ and K . However, in the difference P„ - P . which 
n n n n 

occurs in Eq. (2.30), the odd powers of e will be cancelled. By using 
the results of Appendix B in [7] to Eq. (2.30), one finds that I(x,e) 
has the asymptotic expansion 
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I(x,e) = W(x,e) + V(x,e) = J e 2j [L.+(log e 2 )G.] F(x) (2.31) 

j=l 3 J 

Here, L. and G. are linear operators defined by Eqs. (4.13)-(4.15) in 

J J 

M «v 

[7], with e replaced by e . The functions g k , g k> t.j, t^., h k , h k in 

Eq. (4.14) are defined by Eqs. (B10)-(B17) in Appendix B of [7] 

2 

when e is replaced by e . 

We now return to determine the coefficients a and 8 in the 

n n 

expansions of a(e) and 8(e), respectively. From Eq. (2.24a) it follows 
that fio(x) is analytic on 0 < x < 1, because ^ is analytic and S(x) is 
analytic with S(0) * 0 » S(l) and S ' (0) * 0 * S* ( 1) . It then follows 
that from the recursive nature of Eqs. (2.24a)-(2.24d), the function 
f nm (x) will be analytic for 0 < x < 1 if both GjF(x) and LjF (x) are 
analytic, provided that F(x) is analytic. This will be true if all of 
these functions g k , g k h k , h k are analytic. The functions g k (x) 
and g k (x) are defined in Appendix B of [7] as 

g(x,e) ” l g k (x) e Zk = {[x-a(e)] + £ 2 S(x)} 1/2 (2.32a) 

k*0 " 

g(x,e) - j g~(x)e 2l! - (Cx-B(e)] + c 2 S(x)) 1/2 (2.32b) 

k-0 * 

where o(e) and 8(e) are defined in Eq. (2.18a) and (2.18b) respectively. 

depends on the coefficients <* k in the expansion 
of a(e). They are singular at x=*0 except for certain values of c^, so 
the appropriate « k should be chosen such that g k (x) can be regular at 
x=0. To relate gj ( (x) with © k , one can square both sides of Eq. (2.32a) 
and equate the coefficients of the same power of e. Similarly, g k (x) 
can be related to 8 k at x*l. Detailed derivation of a k and 8 k is 


The functions g k (x) 



available in Appendix A. Using the results Eqs. (A9) and (A17), the 
expressions for a(e) and 8(e) are given by 
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«(e) ■ 4~ 


e - 


16 


(c x c 3 + 6 8 

e + e + 0(e ) 


C 1 2 c l c 2 4 VU 1 u 3 T “ , r*2 / 6 


64 


x _ , d l 2 , d l d 2 4 (d l d 3 + 2d l d 2 } 6 , ft( 8, 
8(e) * 1 - -j=- e + -jg- e p e + 0(e ) 


(2.33) 

(2.34) 


The functions h k (x) are defined in terms of the square root of g k (x). 

Once the a k are chosen appropriately to make g^ (x) analytic, then it 

implies that h k (x) are also regular. Similarly, the h k (x) are also 

regular at x*l, since 9 k (x) are made analytic. Therefore, LjF(x), Gj(x) 

and f nm (x) are analytic provided that F(x) Is analytic on 0 < x < 1. 

2 

The asymptotic expansion of ♦(x,r ;e) can be obtained by using the 
appropriate f nm (x), a(e) and 6(e) which are determined by Eqs. (2.24a)- 
(2.24d), (2.33) and (2.34) respectively in Eq. (2.25). Once the 
potential of the flowfield is found, the aerodynamic properties of any 
given axi symmetric slender body can be evaluated. 


2.1.2 Line Distribution of Singularities Inside a Two-dimensional 
Airfoil 

In Sec. 2.1.1, the potential flow past an axially symmetric slender 
body has been described and the same technique can be applied to the 
two-dimensional airfoil. This idea was presented earlier in [10,11]. 

In this section, the problem will be analyzed in terms of the functions 
of a complex variable. Similar assumptions are made here as in Sec. 
2.1.1. Consider the two-dimensional potential motion of an ideal fluid 
past an arbitrary airfoil which is shown in Fig. 2.2. The problem is 



0 


1.0 



L : body length 


Fig. 2.2 Uniform stream past a two-dimensional airfoil with 
an arbitrary angle of attack. 
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analyzed in Cartesian coordinate (x,y) with the origin at the leading 
edge of the airfoil and the trailing edge is some distance intercepts on 
the x-axis. Now, one should recall the relationship between the 
potential, the stream function and the velocity components in two- 
dimensional motion. As before, the potential is <t> and the stream 
function is i|>. The combination of <j> and into an analytic function as 
a complex potential <|>(z), is given by 

4>(z) - 4>(x,y) + i*(x,y) (2.35) 

and the x-y plane refers to the z-plane (z = x + iy). 

For this study, it is convenient to think of $ as either the real 
or imaginary part of a function 

♦( 2 ) ■ *°(z) + * b (z) (2.36) 

which is analytic in the z = x + iy plane outside the profile curve G of 
the airfoil. The potential Is the potential of the incident flow on 
the airfoil, while «|> b is due to the presence of the body. The 
potential <j> b is represented by the superposition of complex potentials 
of point sources and vortices distributed along an appropriate arc 
inside the airfoil. The boundary conditions of this problem are 

no penetration condition: v<j» . n » 0 on G 

infinity condition: <|> b * tan -1 (y/x) + 0(l/r) as r infinity 
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where n is the unit outward normal vector to the surface and r is the 
total circulation about the airfoil. 

The equation of the profile G of the thin airfoil is given by 
y * e[C(x)±/S(x)] on 0 < x < 1, and e is defined as slenderness 
ratio. This is assumed to be small for a thin airfoil. The leading and 
trailing edges intercept x-axis at x=0 and x*l respectively, and the 
max.|C(x)±S(x)| = 1. It is assumed that both C(x) and S(x) are analytic 
on 0 < x < 1, with C(x) and S(x) vanishing at x=0 and x=l. These 
functions can be expanded in a Taylor series at the endpoints as 

C(x) = l c x n , C(x) = l c (l-x) n 
n=l n n=l n 


S(x) - l 
n=l 

t 0^(0) 
where c = ■■ ■ ■■■ , v / 
n n! 





. S(x) = l d ( X— x ) n 
n-1 " 

n! 


(2.37) 


:(n) 


nT 


ioi 


n! 


/v 

with dj * 0 * dj. 


(2.38) 


In order to solve for it is convenient to find a function, Eq. 
(2.35), which is analytic in z * x + iy outside G and satisfies the 
condition 


Re {[i- e (C'±S’)/2 /S]||} = 0 


on G 


(2.39) 
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Equation (2.39) is verified and explained in Appendix B. Here, <j> b 
represents the superposition of the potential due to point sources and 
vortices distributed along an arc inside the slender airfoil. Thus, 

h 1 p (e) 

* b (z) 3 ’ 27 / ln(z-s) fU.e) d? (2.40) 

a(e) 

where f(s,e) is the unknown source strength density, which is a complex 
expression. The constants o(e) and 8 (e) determine the extent of the 
source distribution and they can be determined after f( 5 ,e) is found and 
they must be inside G. 

Using Eqs. (2.39) and (2.40) in Eq. (2.36), one obtains 
Re {[i -e(C'±S'/2/S)] C X+ MC±/S)] } 

iw 

■ Re (Jj [1-e(C , ±S'/2/S)] / > (*•«) 

Equation (2.41) is a pair of equations, corresponding to where the upper 
(plus sign) and lower (minus sign) surfaces of G, is obtained. By 
adding and subtracting these two equations, a pair of coupled linear 
integral equations can be obtained. Upon denoting x+ie(C+/S) and 
x+ie(C-/S) by xpu and xpl respectively, one has 

d+°(xpl) 
dz - 


2, Re (1(1+1 rf')C £*1< X P U > + ^P 1 > & d»°(xpu) 

dz dz J “ 2/S L dz 

= Re (i 2(x+i eC-g)(l+i eC 1 ) + e 2 S' 
a (x+ieC-£) 2 + e 2 S 


f(5,e) m 


(2.42) 



and 

2, Re {1(l + 1sC')/S[£|-( x P u >-^< x ' )1 >] - f S'[g!^ xpu ' + i|° (x P 1) ]} 


Re t . j 6 2S(l+1sC) - S'(xtieC-s) ~ U)£) d?) 
a (x+ieC-C) + eS 


(2.43) 


In Eq. (2.42) and (2.43), f(5,e:), a(e) and 0(e) need to be determined. 
Both o(e) and 0(e) are found to be power series in e of the form [11] 


o(e) 



6(e) =1-1 $„e" (2.44) 

n*l " 


Using the results of [10] as an instructive example, the function 
f(x,e) shall be of the form 


f(x.e) = . f ( x ’ e > (2.45) 

/C6(e) - x][x - «(c)] 

where f(x,e) is singular in x at x 5 <*(e) and x * 0(e). The function 
f(x,e) is be assumed to be analytic before finding its asymptotic 
expansion about e - 0. Following the procedures given in [10,11], both 
sides of Eqs. (2.42) and (2.43) can be expanded with respect to e, 
without taking account of the dependency of f on e. Since <t>° is an 
analytic function of x + i e[C(x)±/S(x)], the left-hand sides of Eqs. 
(2.42) and (2.43) can be expanded in a power series in e. Then, Eqs. 
(2.42) and (2.43) become 
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00 00 

l e j K.(x) = Re { I e j Lt [f(x,e)]} (2.47) 

j=0 J j=0 J 


Here, Jj(x) and Kj(x) are defined as 


J.(x) ■ 2 ir Re < 
J 1 

cj- Zk (x)s k (x) 

{ J k=0 a 



C(J-i)/23 , , , , 2k k 

2C‘(x) I ( J ,‘) C J 1 Zk (x)S k (x) 

k=o " 


+ S'(x) 

C(j~2)/2] . . i 2k 1 k "1 ) 

J 0 <2kll> CJ 2k - 1 (x)S k (x)J j 

(2.48a) 

K . (x) * 2 tt Re I 
J i 

!«♦,(*) C( 3 T^ +1 ) cJ- 2k - 1 (x)S k+1 (x) 

( J k=0 ^ 1 


X 

H 

1 

1 

r C(j-2)/2] , , , n , k+ , 

1 2C* (x) l (0-l) C J- a '- z ( x )S lt+1 (x) 

L k=0 


+ S'(x) 

C(J ‘j )/ %V'- 1 - 2k (x)S k (x)lj 

k=0 

(2.48b) 

where 4» - (x) » 

J 

( 3T ) {( Hz> *° (z)} z = x 

(2.48c) 


In Eqs. (2.48a) and (2.48b), the notation (?) = n!/(j!)(n-j)! has been 

J 

used and [n] is the greatest Integer not exceeding n. Applying the same 
idea of expanding integral operators asymptotically as in Sec. 2.1.1, 
one can expand the right sides of (2.42) and (2.43). To do this, 
consider the integral operators I p (x,e), p= 0,1 applied to a function 
F(?)/[(0-?)(c-<*)] 1 ^ and F($) is independent of e. The integral 
operators i p are defined as 


I °( x , e ).j B 2[x+i eC(x)-c][l+i eC 1 (x) ] ♦ e 2 s'(x)FU) d? 
a {[x+ieC(x)-c]^+e 2 S(x)> /( 3 - c)(t - a) 


(2.49a) 
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iVc)- , l S ^S(x)[U1eC(x)] -S'(x)Cx.UC(x)- ; ] F( t ) dj (2>4Sb) 
a {[x+i eC(x)-c]^ + e S(x)} /( 0 - ?)($ - o) 

The constants a(e) and 3(e) are defined in Eq. (2.44), where the 
coefficients a n and $ n need to be determined. The operator I p (x,e) can 
be expanded by following the similar method as given in [7] and by using 
the results derived in [11], one obtains 

I P (x,e) = l S L p [F(x)] (2.50) 

q=0 q 

The linear operators L p , p a 0,1 are defined in Appendix C. From the 
definitions given in Appendix C, one has two simple operators 


ir 


CF(x)] = -4 ; (s1n 2 e - X )' 1 [F(s1n 2 9) - F(x» de 

j 0 

(2.51) 

Lj [F(x)]-2, /^F(x) 

(2.52) 

To solve Eqs. (2.46) and (2.47), an asymptotic solution 

for f(x,e) is 

assumed as 


f(x,e) ‘ l f„(x)e n 
n=0 n 

(2.53) 

Here, the functions f n (x) are to be determined. By substituting Eq. 

(2.53) into Eqs. (2.46) and (2.47), one obtains 


I e J J.(x) =■ Re [l j j ^ L? [Mx)] 

j=0 J [_ j=0 n*0 J n J 

(2.54) 
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I Mx) =■ Re I j 0 J+n Li Cf n (x)] (2.55) 

j=o 1 y*o n=o a " j 

By equating coefficients of the like powers of e in the corresponding 
equations, Eqs. (2.54) and (2.55), finds 

(2.56) 

(2.57) 

(2.58) 

(2.59) 

Applying the results of Inversion of operators In [10] with Eqs. (2.58) 
and (2.59) in Eqs. (2.56) and (2.57) respectively, one has 

Im f (x) = f (x) + Ira f (0) (2.60) 

n n n 

where 

V*> -i{ + L 3° [f n-j (x)] l 


I™ Uj Cf n (x)]> - -0 n (x) - Lj° tf n . J (x)]J 

Re <l£ [f n <*)]> - K n (x) - Re l] [f n _j<x)]J 
Using the results of Eqs. (2. 51) and (2.52), one finds that 

Im {Lg[f n (x)]} = Im { L° [f n (x) ]} = L° [ Ira f p (x)] 

1 Im [f (sin 2 e)] - Ira [f_(x)] 

= -4/ ^ H de 

0 (sin^e - x) 

Re {l£ cyx)]} - 2, J Re [f„(x)3 
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IT 

I - -- f - |J„(S) + I" 

* 0 (sirre - x) L n j 


[ n n 15 = si n 2 e 

J (5) + Im l LTCf n ,•(€)] d0 (2.61) 

j=l J J J E = X 


and 


Re f »w- i? { V*> - «• Ij L ] J » 


> 0 (2.62) 


S=a 


In Eq. (2.61), the notation [G(s)] * G(a)-G(b). From Eqs. (2.61) and 

(2.62), the function f p (x) can be determined recursively, once the 
constants Im f n (0) are found. In order to determined Im f n ( 0 ) , one has 
to use the condition that 


where 


8 ** 

Im / f( 5 ,e) d? = r 

a 


r - I r. e J 


(2.63) 


j=0 


Here, each i\ is a prescribed constant Independent of e. By using Eqs, 
(2.45), (2.53) and Eq. (2.61) In Eq. (2.63) can be expressed as 


f i (5) etc • . 

„ - , *—■ - — + » I Im fj(0) 

j’O <1 /(B - c)(c - o) J"0 3 


• * 3 

r - I J I 


(2.64) 


By setting 5 * (a-a) sln^e + a and expanding each integrand in Eq. 
(2.64) in Taylor series about e ■ 0, one obtains 

r - 2 l E j Im l ±.f {(|-) i f . [(a - o) s1n 2 e + a]} da 
j=0 n=0 n ' 0 36 J ' n e=0 


+ ir l Im f. 
j=0 J 


(0) 


(2.65) 
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Equating coefficients of like powers of e on each side of Eq. (2.63) and 
(2.65), yields 

IT 

Ini yo)--^- | J L- /{(f-) k f n _ k C(e-o.)s 1 n 2 0 + «]) de (2.66) 

n ir * k=0 0 de n k e=0 

Since f Q (x) depends only on Jg(x), the Im f n ( 0) can be determined 
recursively from Eq. (2.66). Thus, functions f n (x) can be found through 
Eqs. (2.61), (2.62) and (2.66) in a recursive manner. 

Determination of the constants a(e) and 8(e) is discussed now. 

From Eqs. (2.61) and (2.62), one knows that Re f g(x) and Im fg(x) depend 
on <j>° (z) which is analytic. Thus, it follows from the recursive nature 
of Eqs. (2.61), (2.62) and (2.66) that all the Re fg(x) and Im fg(x) 
must be analytic for 0 < x < 1 if the Lq[F(x)] is analytic. This will 
be true if and only if each gj(x) and hj(x) in Eqs. (C2) and (C3) is 
analytic on the interval 0 < x < 1. In the manner of Appendix A, one 
can determine a k uniquely by demanding that g(x,e) is regular at x=0. 
Similarly, all of the 8 k can be determined by requiring that the 
function h(x,e) is regular at x»l. By squaring both sides of Eqs. (C2) 
and (C3) and equating the coefficients of the same power in e, one can 
find out all the appropriate and e k , such that singularities at x=0 
and x=l of functions gj(x) and hj(x) are eliminated. The similar 
technique for the determination of a(e) and 8(e) is given in Appendix A, 
and the general expressions for a s< and 8 k are available in Appendix C. 
Using the results In Appendix C and the definitions of Eqs. (2.37) and 
(2.38), the leading terms for o(e) and 8(e) are given by 


, . d l 2 i Cld 
«(e) ■ 4- e - 


1 U 1 3 d l (d 2 + 4 . n , 5 X 

— e 16 e + 0(e ) 


16 


(2.67) 
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„/ % , d l 2 1 c l d l 3 t d l (d 2 + 40 1 } 4 , n , 5, 

8(e) * 1 - -f- e - e + jg e + 0(e ) 


( 2 . 68 ) 


Once f(x,e) has been obtained, the expression for 4 >(z ) is obtained by 
using Eqs. (2.45), (2.53) in (2.40) and (2.36) is given by 


4>( z »e) = $°(z) 


1 

• "27 


r» K 

I e 

k=0 


8(e) In(z-s) f k (c) d 5 
•(e) / CS(e) - ?][? - «(•)] 


(2.69) 


Equation (2.69) is the desired uniform asymptotic expansion for <j>. 

Suppose now that the thin airfoil has a sharp trailing edge such 
that near x=l, S(x) has the expansion 

S(x) ■ l d.(l-x) j , d 2n+1 # 0 n > 1 (2.70) 

j*2n+l J 1 

where d. = (-1) J S^ J ^(l)/j. From Eq. (2.70), it follows that one can 

J 

take 8(e)=l, the asymptotic expansion obtained above remains valid when 
Eq. (2.70) holds. The mathematical proof is available in Sec. 7 of 
[10,11]. Physically, the velocity of the fluid near the sharp trailing 
edge should be finite. As the velocity of the fluid with the potential 
given by Eq. (2.69) near the trailing edge Is examined, the first term 
in Eq. (2.69), 4 >°(z) obviously gives a finite velocity at the trailing 
edge, because <j»°(z) is analytic at z»l. The remaining terms in Eq. 
(2.69) contribute velocity components of the form 


1 



1 


z - S 


F(e) 

/[C - «(e)] (1-7) 


d? 


(2.70a) 


Equation (2.70a) remains bounded as z>l, if and only if F(1)=0. 
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Therefore, one shall make each f k (x), k > 0 vanishes at x=l by properly 
choosing the arbitrary constants of Im f k (0). The mathematical proof of 
choosing Im f k (0), so that each f k (x), k > 0 vanishes at x=l, Is shown 
in Sec. 7 of [10], i.e., 

Im f k+l (0) “ “ f k+l (1) (2 ’ 71) 

If the flow described by $° is a uniform stream, the circulation 
can be determined by Eq. (2.65), such that the lift L(e) of the airfoil 
is given by 

L(e) = p U r(e) (2.72) 

00 

where p is the density of the fluid and U is the speed of the uniform 

oo 

stream at Infinity. By using Eqs. (2.65) and (2.71), one finds that Eq. 
(2.72) becomes 


IT 



-irf.(l) (2.73) 

The uniform asymptotic expansion for <|>(z ; e) can be determined by finding 
f n (x), a(e) , and 3(e). Then, the velocity, pressure distributions on 
the surface, and also the lift coefficient of an arbitrary airfoil can 


be determined. 
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2.1.3 Second Order Joukowski Airfoil 

Certain assumptions have been made in analyzing the problem in Sec. 

2.1.2. The upper and lower surface of the thin airfoil have the 

equations y = e[C(x)+/S(x)] and y= e[C(x)-/S(x)] for 0 < x < 1; with 

C(0)=S(0)=C(1)=S(1)=0. The equations of the thickness and the camber 

line of the airfoil are y 3 e/S(x) and y = eC(x), respectively. 

Here e is the slenderness ratio, while S(x) and C(x) give the thickness 

and camber distributions. Varying e produces a family of affinely 

related profiles. More examples are given by Van Dyke in [1]. Now one 

expects to seek an expression for a thin Joukowski airfoil by using the 

perturbation method. The parametric form of the Joukowski airfoil is 

easily obtained by means of conformal mapping [15]. One employs a 

transformation x(s,n), y(£,n), where the s-n plane refers to as 

the c -plane which will map the region and its boundaries in 

the 5 -plane into another desired configuration in the x-y plane, which 

is referred to as the z-plane. Fig. 2.3 shows the notation of the 

Joukowski transformation of a circle into an airfoil. The center of the 

i 6 

circle with radius a at any arbitrary point is given by y = m e . The 

-i a 

symbol y Q denotes the intersection of the vector a e p with the y- 
axis. A circle with center at y Q and radius of magnitude |y Q - c-j-| 
transforms Into a circular arc, which is considered as the skeleton of 
the airfoil. The arbitrary circle, the points of which are described by 

ie 

C s a e + y (2.74) 

transforms into an airfoil. The angle 8 determines the mean curvature 



o 


2C 


Fig. 2.3 Illustrating the derivation of an arbitr 
airfoil . 
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of the profile and the magnitude |u - u 0 | determines the thickness of 
the profile. 

The circle is chosen such that it passes through ? = Zj * C, where 
C is real positive number. The corresponding point on the airfoil is 
given by 

C 2 

Z(e)=C+£- (2.75) 

It should recall that the center p of the circle. Its radius a, and the 
angle 0 are connected by 

-ie 16 

C * a e + m e (2.76) 

Since C is a real positive number, one can separate the real and 
Imaginary parts of Eq. (2.76) as follows: 


C = a cos 0 + m cos5 

(2.77a) 

a sins = m sin5 

(2.77b) 


Using Eq. (2.74) in Eq. (2.75), with the expressions given by Eqs. 
(2.77a) and (2.77b), Z(?) can be separated into real and imaginary parts 
as x and y respectively. After further simplification, the parametric 
form of the airfoil is obtained as follows: 

x = a(cose + h) I 1 + -5 5 i 1 (2.78a) 

L a^(l + h* + k* + 2 h cose + 2 k sine)J 
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y = a(sine + k) 1 - 9 1 (2.78b) 

j_ a ^(1 + h‘ + k* + 2h cose + 2k sine)J 

where h = (m/a) cosfi , k * (m/a) sin6 (2.79) 

The tral ling-edge point of the airfoil is denoted by Zj with e * -e in 
Eq. (2.75), one finds x = Re(Z-j-) * 2C and y = Im(Zy) = 0. Thus, the 
trailing edge of the airfoil intercepts the x-axis at the point T with 
0T=2C. In order to have an asymptotic expression for a thin Joukowski 
airfoil a small parametric, should be introduced to expand the 
available equations in power series. Generally the quantities h, 
k, 0, and m should be small. If we choose m=e.jC, then Eqs. (2.77a), 
(2.77b) and (2.79) can be expanded to the second order of accuracy as 
follows: 

a “ cosj ( X " e i cos5 ) (2.80a) 

h » (cos 0 cosfi) + 0(6^) (2.80b) 

k * (cos0 sinfi) e-j + O(e^) (2.80c) 

Using Eqs. (2.80a)-(2.80c) in Eqs. (2.78a) and (2.78b), and retaining 
the terms with order of e-j , one obtains 

C 2 2 

x = {cosed + cos 0 ) + ejt-cossdose + cos0)(l + cos 0 ) 

o 

+ 2 cos 0 cose(cosfi - cos 0 cosfi cose - cose sinfi sine)]} 
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+ 0(ej ) 

(2.81a) 

2 1 

y * Csine cose ej[2cos e cose - (cose + ^j)] 


2 

+ Csine + Ce 1 sin5[sin 2 e + 2sin 2 e cos 2 e] + 0(e 2 ) 

(2.81b) 

o 

By examining carefully the order of magnitude of sine, cose, sin e 

2 2 
and cos e, one finds that both cose and cos e are 0(1) and sine and 

2 2 

sin e are 0(e 1 ) and O(e^) respectively. Neglecting the terms of 
0^), Eq. (2.81a) can be further simplified as 

x = 2Ccose + O(e^) 

(2.82) 

Furthermore, Eq. (2.81b) can be simplified by comparing the 

2 2 

magnitude of sine, sin e, cose and cos e as mentioned. To 
of 0(e^, Eq. (2.81b) becomes 

order of 
retain terms 

2 2 

y * 2CejCos6 sine (cose - 1) + 2Ce^sin6 sin 8 + O(e^) 

(2.83) 

Introducing the nondimensional quantities x = x/2C, and y = 
(2.82) and (2.83) become 

y/2C, Eqs . 

x * cose + O(e^) 

(2.84a) 

o p 

y = GjCosS sine(cose - 1) + e^sinfi sin 8 + O(e^) 

(2.84b) 


Upon eliminating sine and cose in terms of x, Eq. (2.84b) becomes 
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y = e 1 (-cos5)(l - x) / 1 - x 2 + sin5(l - x 2 ) 

+ 0(s 2 ) -1 < x < 1 (2.85) 

Obviously, if 6 » ir, Eq. (2.85) will yield a symmetrical second order 
Joukowski airfoil, while the second term vanishes. It is consistent to 
the second order Joukowski airfoil shown in Sec. 4.6 of [1]. It also 
implies that the second term represents the equation of the camber 
line. Therefore, Eq. (2.85) can be rewritten as 

y = y s + y c (2.86) 

where y g * e 1 (-cos6)(l - x) / 1 - x 2 and y c * sinfi(l - x 2 ) . 

To make Eq. (2.86) consistent with the required form of equation, simple 
translation and scaling of the coordinates are done, such that Eq. 

(2.86) becomes 


2(-COS8) / x(l - x) 3 


(2.87b) 

2 sin5 x(l - x) 

0 < x < 1 

(2.87b) 


By defining e * 2(-cos6)ep with e Is the slenderness ratio, Eq. (2.87b) 
yields y ■ e /S(x) and y « e C(x) 

w U 

where 

S(x) » x(l - x) 3 


(2.88a) 
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C(x) * Ax(l - x) , A = -tanfi (2.88b) 

S(x) and C(x) give the thickness and camber distributions of an airfoil 
respectively, and A is related to the airfoil camber. 

2.2 Other Solution Techniques 
2.2.1 Some Exact Solutions 

The analytical solutions of potential flow problems associated with 
the motion of a solid body through an ideal fluid are considered here. 
The final aim is to determine the pressure distribution over the surface 
of the body and the lift force acting on the body whenever available. 

The analytical solution for elliposidal body of revolution and two- 
dimensional elliptic airfoil are obtained in this section. 


1. Analytical Solution for Ellipsoidal Body in Uniform Stream 

Consider a rigid ellipsoidal body moving through a uniform ideal 

fluid. The flow is assumed to be steady, incompressible and irrota- 

tional. To simplify the problem, a body -fixed reference frame and 

elliptical coordinate system U.n.B) are introduced. The derivation is 

described in Appendix D. Using the results of Eq. (D5a), one is able to 

obtain 0 = U c e c + U e where U = rX- -§* and U = - tX- || . 

5 5 n n 5 h h an n h h as 

n P S P 

Here, U_ and U are velocity components in £ and n directions 
% n 

respectively, and 


2 2 

a i|» a i|> aip a^ _ n 

^* 17 ' coth5 « ’ cot '’ " 


(2.89) 


Consider the Neumann boundary condition on 5 ® and the infinity 
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condition as £ tend to infinity. Making reference to Appendix D, the 
boundary conditions can be simplified as follows: 


♦(5.n) = C' on £ ■ £ r 


(2.90a) 


u A ? 9 

i|>(£,n) + — sinh £ sin n + C' ' as £ + infinity (2.90b) 

Equation (2.89) along with Eqs. (2.90a)-(2.90b) can be solved by 
separation of variables. One finds out that the velocity component 
in e^ direction on the body surface is zero (due to the no penetration 
on the surface of the body), and the velocity in e^ direction is given 
by 


v - 


U si nn 

oo 

coshe - 1 


[ o cosh £ - 11 / y 5“ 

cosh£ o + sinh c £ o ln(- - sin j - - ) /cosiT£ 0 - cos n 


(2.91) 


With application of Bernoulli equation, the pressure distribution on the 
surface of the elliposidal body is given by 


U s 2 

c d - 1 - ( TT> 


(2.92) 


2. Analytical Solution for Elliptic Airfoil in Uniform Stream 

In Sec. 2.1.3 the idea of conformal mapping has been introduced to 
get the shape of an airfoil. Here, the pressure distribution over the 
airfoil, by means of Joukowski transformation is discussed. This method 
can be used directly to the actual problem of by means of mapping the 
flow past an arbitrary circle, for which the solution is already 
known. To define the arbitrary circle, we describe the points by 
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t 


(2.93) 


where r Q = (a+b)/2 with a and b are the major and minor semi -axis 
lengths of the ellipse respectively. Making reference to [15], the 
complex velocity on the airfoil surface is given by 

-ie 

i 2U sine e 

W( Z ) =_ (2.94) 

where dz/d 5 * 1 - C 2 /^ 2 and C = (v4 2 - b 2 )/2. After simplification, 
the velocity distribution along the contour of the ellipse is given by 


u s 00 
u 

00 


2 2 1/2 

(l+k)(a -xj 

2 2 2 T/2 

[a^ + x^(k -l)] 1 ' ^ 


(2.95) 


where k = b/a. The pressure distribution on the elliptic airfoil can be 
computed through the Bernoulli equation. 


2.2.2 Panel Method 

A considerable number of panel methods has been developed in the 
past two decades to analyze the steady, inviscid and irrotational flow 
fields past a body of arbitrary geometry. The Incompressible flowfield 
is governed by the Laplace equation which can be transformed to an 
Integral equation relating the perturbation potential to the source and 
doublet singularities distributed over the surface. There are different 
discretization techniques among panel methods for solving the Integral 
equation. Basically, all the methods approximate the surface by 
elemental panels of prescribed geometric shape and singularity variation 
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and the Integral equation is solved by enforcing a boundary condition at 
a control point or points on the panel. There is a summarized table in 
[16], showing several useful panel methods with panel geometry, 
singularity variation and boundary condition specification. In this 
section the governing equation and the basic idea of discretization are 
revi ewed . 

Consider steady, inviscid, irrotational and incompressible fluid 
flow in a bounded domain D, Fig. 2.4, governed by the Laplace equation 

V 2 «j> = 0 (2.96) 

where 9 * 0^ + v<f>. The perturbation potential <fr(P) at any point P can 
be expressed as the potential induced by a combination of source (a), 
and doublet (u) singularities distribution over a bounded surface S 

♦<p> - |/ .<- ^r) ds ♦ „ Is- <£) dS (2.97) 

where r is the distance from the point P to the bounded surface 
and n represents the normal to the surface S directed into the fluid 
domain. The boundary condition in this study is the flow tangency 
condition or 

9 . n = 0 (2.98) 

Since there are two unknowns a and y, while there is one boundary 
condition, there exists an infinite number of source and doublet 
combinations which satisfy Eq. (2.97) and one can specify the doublet 
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strength and solve for the source strength or vice-versa. Using the 
idea of real external and imaginary internal fields introduced in [17], 
the condition for the perturbation potential interior to the region is 
constant and set equal to zero, i.e., 

= 0 or ♦ * u (2.99) 

Thus, it follows that 

a * V<fr . n (2.100) 

For Neumann problems, the source density distribution is given by 
Eq. (2.100), the doublet distribution is to be determined. For 
Dirichlet problems, the doublet distribution is known from the 
prescribed potential on the boundaries, and the source distribution is 
to be determined. 

In order to Introduce lift, wake networks extending from the 
trailing edge to downstream infinity are considered. As the free stream 
Reynolds number increases, the thickness of the wake region diminishes 
and is often assumed to be zero. Since fluid properties are continuous 
across the wake, the wake region is usually modelled as a stream surface 
having zero pressure loading Imposed by the surrounding flow. The wakes 
can be represented by a doublet sheet with zero thickness and streamwise 
doublet gradient, emanating from the trailing edge and parallel to the 
free stream. Determination of the strength of the wake constitutes the 
Kutta condition which is available in [16]. Some typical numerical 
schemes were already presented in the literature [16-21]. 
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Chapter 3 

APPLICATION OF PERTURBATION METHOD 

3.1 Perturbation Technique for Slender Bodies of Revolution 
The formulations of Sec. 2.1.1 are applied to the case of a uniform 
ideal fluid flow past an axi symmetric slender body. The coordinate 
system is fixed to the body and the velocity of the stream at infinity 
is U^. The Incident stream function is given by 

*°(>c,r 2 ) = U r 2 /2 (3.1) 

QB 

The stream function y° can be expanded in powers of e 2 and the 
coefficients i|>.(x) can be obtained by using Eq. (3.1) in Eq. (2.19b); 

J 

this yields 


♦ X (x) = UJ2 , *j(x) =0 for j > 2 (3.2) 

By inserting these coefficients into Eqs. (2.24a) - (2.24d) together 
with the definition of operators Lj and Gj, the first six leading terms 
of f nm (x) are given as follows: 

f l 0 (x) - » S’ 

f 20 ( x ) s ” 7 3x (L l f 10^ 


(3.3a) 
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’-i u .37ff--ilr- ss “ +s “ s 

1-x 0 

+ S / [S'(x+v) - S'(x) - vS"(x)] v" z dv 

0 

- s ; X [S'(x-v) - S' (x) + vS"(x)] v“ 2 dv} (3.3b) 

0 

f 21 (x) “ " 7 If (G l f 10 ) " I U « Hx {SSB) (3,3c) 

f 30 (x) = 0 (3.3d) 

f 3 l( x ) “ - 757 ( L i F 2l + G 2 f 10 + G 2 f 2(P 

f 2 0 ) -Ur( $2 f i6‘) 


- ST Cfh< S log [i^] - S> + (f - ^f*] 


+ S / [f 21 (x-v) - f 21 (x) + vf' l( x)] v~ 2 dv 


1-x , 


- s J Cf 21 (x + v) - f 21 (x) - vfjj (x)] v" z dv} 

(3.3e) 

f 32 (x) = < S l f 21 ) = U » H7 [S ( SS " )] 

(3.3f ) 


Equations (3.3a) - (3.3f), together with (2.18a) and (2.18b) depend upon 
the profile curve. The expressions for f nm (x), a(e) and e(s) can be 
evaluated If S(x) Is specified. The complete asymptotic expansion 
for 4>(x,r 2 ;e) can be obtained from Eq. (2.25) with » Ujc. The flow 
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velocity components over the body surface can be determined by taking 
appropriate derivatives of <t>. Next the Bernoulli equation is used and 
the pressure coefficient on the body surface is evaluated as 

U s 2 

Cp * 1 - (/) (3.4) 

00 


where U s is the velocity distribution on the surface of the body. The 
profile curve is defined by r » e /S(x) on 0 < x < 1. 

In this study, two different axi symmetric bodies are analyzed, and 
their profiles are given as follows: 

1. Ellipsoidal Body 
S(x) = 4x(l - x) 

2. Dumbbell Shaped Body 

S(x) * 4bx(l - x)[l - bx(l - x)] , b > 2 

With the profile shape r a e /S(x) on 0 < x < 1 mentioned, one can 
obtain explicit expressions for f nnj (x), a(e) and 0(e) from Eqs. (3.3a)- 
(3.3f), (2.18a) and (2.18b) respectively for each of the above two cases 
(Figs. 3.1 and 3.2). Lists of the required expressions are presented in 
Appendix E. 


3.2 Perturbation Technique for Thin Airfoils 
The formulation given in Sec. 2.1.2 is applied to the flow problem 
when the airfoil is at rest in an ideal fluid which has unit velocity at 
infinity. The airfoil is at an arbitrary angle of attack y to the x- 
axis. Thus, the complex potential for the incident stream is given by 


♦°(i) ■ e' 1 ’' z 


(3.5) 
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Introducing Eq. (3.5) in Eq. (2.48c), one obtains 


Mx) 


if j = 0 
if j >1 


Using Eq. (3.6) in Eqs. (2.48a) and (2.48b) yields 


Jq(x) * 4ir siny 
Jj(x) * -4ir C'(x) COSy 


Jj(x) - 0 


K l(x) * -2ir S'(x) COSy 


K j ( x) * 0 


if j > 2 


if j = 0,2, .., 


The functions f n (x) are obtained by inserting these coefficients 
together with the expressions in Appendix C defining the operator 
Lj, into Eqs. (2.61) and (2.62). The three leading terms f n (x) are 
computed as follows: 


Re f 0 (x) = 0 


f q ( x ) « 2x siny 


!m f Q (x) = — - siny 


Re f 


1^ = " 2ir S(x)^ ^irS 'cosy + Re Lj [f Q (x)]J 


Cf 0 (x)] - 4 tt C'(x) cosy 
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7 2 

X f COS 0 

" T J 9 

x 0 sin 9 - x 


[ Ira L ? 


Cf 0 (?)] - «w C'(5) COST 


rl 


2 

5=sin 0 


J5=x 


d0 


ir 

r y ~ 

Im f , (0) = — - — J f.(sin 2 0) d0 

JL TT TT q Jl 


(3.9c) 


Re f 2 (x) =■ - ^Re lJ [fjfx)] + Re t\ [f 0 (x)]l (3.10a) 


f 2 U) - $7 


7 {*" *-i Cfj 


(x)] + Im 1-2 [f 


o^l 


(3.10b) 


x 

T 

x 


X 

J <**» .., L 1° Cf.( 5 )] + to L° [f 0 (5)]l 
0 slnVx L 1 1 1 “ J 


2 

5=sin 0 
S=x 


d0 


x 

r *2* /w 

Im f 2 (0) “ “ “ 7 J [f 2 ( sin20 ) +2 siny(a 2 cos 2 0 - 0 2 sin 2 0)]d0 (3.10c) 


The operators L^, l|, L®, and L 2 are given in Appendix C. When Eqs. 
(3.8a)-(3.10c) are used in Eqs. (2.53) and (2.45) they give the 
asymptotic expansion for f(x,e) up to 0(e ). Once f(x,e) has been 
found, one can obtain an expression of <fr(z;e) by using Eq. (2.69) and 
the pressure coefficient can be evaluated from Eq. (3.4). Since Eqs. 
(3.8a)-(3.10c), (2.67) and (2.68) are functions of S(x) and C(x), once 
they are specified, the complete expansion for <f»(z;e) up to 0(e ) can be 
determined. 

In this study, elliptic airfoils and second order Joukowski 
airfoils have been analyzed using the above method. For a symmetric 
airfoil the function C(x) is set to zero; only the thickness solution 
has to be solved. If the airfoil is symmetric and is set at zero angle 
of attack, then there is no circulation about the airfoil and all r n are 
set equal to zero in determining the constants Im f n (0). Moreover, if 
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the airfoil has a sharp trailing edge, special treatment is applied to 
obtain each r n> which is given in Eq. (2,71). The profile curve of the 
thin airfoil Is defined by y * e[C(x)±/S(x)] on 0 < x < 1, with small 
slenderness ratio e. Three different airfoils are used as test cases. 
Figs. 3.3, 3.4 and 3.5, and their geometries are as follows: 

1. Elliptic Airfoil 
S(x) = 4x(l - x) 

2. Symmetrical Second Order Joukowski Airfoil 
S(x) ® x(l - x) 3 

3. Cambered Second Order Joukowski Airfoil 
C(x) * Ax(l - x) 

S(x) = x(l - x) 3 

Here, A is related to the airfoil camber and is defined in Eq. 

(2.88b). Using the above Information, expressions for 
f n (x), r n , o(e) and 0(e) can be derived in each case. A brief 
derivation and the lists of the resulting expressions are included in 
Appendix E. 




Fig. 3.3 Elliptic airfoil of 10% thickness. 
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Fig. 3.4 Symmetrical Joukowski airfoil of 6.5% thickness. 
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Fig. 3.5 Cambered Joukowski airfoil of 6.5% thickness. 
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Chapter 4 

RESULTS AND DISCUSSION 

In this chapter an assessment of the merits of the perturbation 
analysis method Is made. The pressure distribution on axi symmetric 
bodies and two-dimensional airfoils obtained by the above method are 
compared with those of panel method as well as exact solutions whenever 
aval i able. The results of the panel method are obtained by using the 
low-order panel method code (program VSAERO, developed by Maskew 
[18]). The analytical expressions for the ellipsoidal body and elliptic 
airfoil are avaliable in the literature presented in Sec. 2.2.1. 

4.1 Axi symmetric Bodies 

Since the analytical expression for an ellipsoidal body is 
available, this body is selected to demonstrate the accuracy of the 
perturbation analysis method with the variations of the slenderness 
ratio and also the number of terms in the expansion of t. Figs. 4.1a- 
4.2c illustrate a comparison of the pressure distributions, as calcu- 
lated by the perturbation analysis method with different order of 
accuracy in the expansion of 4 , with analytical solutions for the two 
values e = 0.1 and e * 0.2, corresponding to L/d = 5 and 2.5, respec- 
tively. The results show that the perturbation solutions agree well 
with the exact solutions to a high order of accuracy, i.e., up to 

Q O 

0(e ). With the same order of accuracy in 4> up to 0(e ), the pressure 
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distributions obtained from the perturbation analysis method and the 
panel method are compared with the exact solutions, with different 
numbers of panels in the chordwise direction for e = 0.10. The 
calculations are presented in Figs. 4.3a-4.3c. It is seen that the 
results using 30x6 panels and 50x6 panels yield better agreement with 
the exact solutions than using 10x6 panels. Again for the same order of 

O 

accuracy In <t> up to 0(e ), and using 50x6 panels, three sets of pressure 
distributions are compared with four different values of e. The results 
are presented in Figs. 4.4a-4.4c. The perturbation solutions agree well 
with the exact solution with small slenderness ratio up to 0.1, and the 
panel solutions compare well with exact solutions for all values of e. 
The dumbbell shaped body is also analyzed by the perturbation analysis 
method with the expansion of $ up to 0(e 4 ), as well as by the panel 
method. A comparsion of pressure distributions obtained by the above 
method and the panel method, using 50x6 panels, with different values 
of e and profile parameter b are presented in Figs. 4.5a-4.6b. For 
small e and b, the perturbation solutions agree well with that generated 
by the panel method. In the present examples, each perturbation case 
required only 2 seconds computing time versus 5 seconds (30x6 panels) or 
14 seconds (50x6 panels) for the panel method. This computation time 
comparison is based on the CDC Cyber 173. 

4.2 Two-Dimensional Airfoils 

The analytical pressure distribution on the surface of the elliptic 
airfoil is avaliable in the literature discussed in Sec. 2.2.1. It is 
convenient to do numerical experiments for the perturbation analysis 
method using this airfoil. Perturbation analysis solutions for an 
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elliptic airfoil, set at zero angle of attack, are compared with exact 
solutions, with different order of accuracy in <j> expansion for e = 0.10 
and e * 0.15, corresponding to the t/C ratio = 20% and 40% respectively, 
as presented in Figs. 4.7a-4.8b. It is apparent that the perturbation 
analysis solution gives quite good agreement with the analytical 
solution for high order of accuracy in $ up to 0(e 4 ). The pressure 
distributions calculated from the perturbation analysis method are 
compared with the exact solutions, with different values of slenderness 
ratio and these are shown in Figs. 4.9a-4.9c. The results show that the 
perturbation solutions agree closely with exact solutions up to 20 
percent thick elliptic airfoil. 

Results of the symmetrical second order Joukowski airfoil obtained 
by the perturbation analysis method with the order of accuracy of 4> 

3 

expansion up to 0(e) are compared with the panel method in the 
following section. Here, 3.25%, 6.5%, 10% and 13% thick airfoils are 
chosen as the test cases. For the calculation presented, the thickness 
of the airfoil, angle of attack, and the number of panels are given for 
each of the calculation. The calculations presented in Figs. 4.10a- 
4.10c and 4.11a-4.11c are obtained using different numbers of panels in 
the chordwise direction for symmetrical Joukowski airfoil of 6.5% 
thickness at 0° and 6° angles of attack, respectively. The results do 
not give a smooth curve near some regions of the leading edge until 50 
panels are used. A comparsion of the perturbation analysis solutions 
for the symmetrical second order Joukowski airfoil at 0°, 6° and 12° 
angles of attack with panel method are presented in Figs. 4.12a-4.15c. 
The perturbation solutions agree with the panel solutions except in some 
region on the leading edge of the airfoil surface. Table 4.1 ill us- 
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trates the appreciable difference between the lift coefficients of the 
two methods from the exact solutions. 

The perturbation analysis method has been applied to the cambered 
second order Joukowski airfoil with the order of accuracy in the 

3 

expansion for $ up to 0(e ). Calculations are presented in Figs. 4.16a- 
4.19c. The results generated by the panel method with 0°, 6° and 12° 
angles of attack are used for comparison. Both pressure distributions 
agree closely for small slenderness ratio and angle of attack up to e = 
0.03248 and y - 6° respectively, except at some regions near the leading 
edge. The lift coefficients obtained from the perturbation analysis and 
the panel method compared with the exact solutions are shown in Table 
4.2 . In the present example, the computation time for each pertur- 
bation case required 8 seconds, while 0.5 seconds (10 panels), 2 seconds 
(30 panels) and 5 seconds (50 panels) were required for each panel - 
method solution. 
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(a) e = 0.1 ; <|> expanded up to e 2 . 


Fig. 4.1 Comparison of pressure distributions from the perturbation 
analysis method with the analytical solution for 
ellipsoidal body. 















ComDarison of pressure distributions from the perturbation 
analysis method with the analytical solution for 
ellipsoidal body. 
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(a) e = 0.1 ; 4> expanded up to e 6 ; 60 panels with cosine spacing. 


Fig. 4.3 Comparison of the perturbation and exact solutions with 
the nanel method for elliosoidal body. 




(b) e = 0.1 ; 0 expanded up to e 6 ; 180 panels with cosine spacing. 


Fig. 4.3 Continued. 
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(a) e = 0.05 ; $ expanded up to e 6 ; 300 panels with cosine spacing. 


Fig. 4.4 Pressure distributions from the perturbation analysis 
and the panel methods are compared with the exact 
solution for ellipsoidal body. 
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(c) e = 0.20 ; <{> expanded up to e 6 ; 300 panels with cosine spacing. 


Fig. 4.4 Continued. 
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(d) e = 0.25 ; <j> exDanded ud to e 6 ; 300 panels with cosine spacing. 


Fig. 4.4 Concluded. 
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(a) e = 0.03 ; b = 2 ; <p expanded up to e*; 

300 panels with cosine spacing. 

Fig. 4.5 Comparison of pressure distributions from the Derturbation 
analysis method with the panel method for dumbbell shaped 
body. 
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(a) e = 0.03 ; b = 3 ; <j> expanded up to e 4 ; 

300 panels with cosine spacing. 

Fig. 4.6 Comparison of pressure distributions from the perturbation 
analysis method with the panel method for dumbbell shaped 
body. 






YC XJ = .1CC»5CRT(4X:1-XI) 


3 



(a) e = 0.1 ; <j> expanded up to e 2 , t/C = 20% . 

Fig. 4.7 Comparison of pressure distributions from the perturbation 
analysis method with the exact solution for elliptic 
airfoil . 
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Fig. 4.8 Comparison of pressure distributions from the perturbation 
analysis method with the exact solution for elliptic 
airfoil . 






(a) e = 0,05 ; <j> expanded up to e 3 ; t/C = 10% . 


Fig. 4.9 Comparison of pressure distributions from the perturbation 
analysis method with the exact solution for elliptic 
airfoil. 











(a) e = 0.03248 ; y = 0° ; <j> expanded up to e 2 ; t/C = 6.5% ; 

10 panels with cosine spacing. 

Fig. 4.10 Comparison of pressure distributions from the perturbation 
analysis method with the panel method for symmetrical 
Joukowski airfoil. 
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X 


(c) e = 0.03248 ; Y = 0°; (j> expanded ud to e 2 ; t/C = 6.5% ; 
50 panels with cosine spacing. 


Fig. 4.10 Concluded. 
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(a) e = 0.03248 ; y = 6°; <J> expanded up to e 2 ; t/C = 6.5% ; 

10 panels with cosine spacing. 

Fig. 4.11 Comparison of pressure distributions from the perturbation 
analysis method with the panel method for symmetrical 
Joukowski airfoil. 



87 



X 

(b) e = 0.03248 ; y = 6 ; <p expanded up to e 2 ; t/C = 6.5% ; 
30 panels with cosine spacing. 


Fig. 4.11 Continued. 









(a) e = 0.01624 ; y = 0° ; $ expanded up to e 2 ; t/C = 3.26% ; 

50 panels with cosine spacing. 

Fig. 4.12 Comparison of pressure distributions from the perturbation 
analysis method with the panel method for symmetrical 
Joukowski airfoil. 
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X 


(b) e = 0.01624 ; Y = 6° ; <p expanded up to e 2 ; t/C = 3.25% ; 
50 panels with cosine spacing. 


Fig. 4.12 Continued. 
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(a) e = 0.03248 ; y = 0° ; <j> expanded up to e 2 ; t/C =6.5% ; 

bO panels with cosine spacing. 

Fig. 4.13 Comparison of pressure distributions from the perturbation 
analysis method with the panel method for symmetrical 
Joukowski airfoil. 
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(c) e = 0.03248 ; Y = 12° ; $ expanded up to e 2 ; t/C = 6.5% ; 

50 panels with cosine spacing. 


Fig. 4.13 Concluded. 


YC X) - .150s! CC X) »■/- SORT! SI XI 11 NrHMr 50 
CC X) = -OCOXI L -XI SI XI =X( l -X) *»3 CN= 50 


T/C= 3.75X ANGLE OF ATTACK:: G.00 CSGREE 


ALFHAz .572CE-G2. 0. 


EETPr! .0000 


CROSS PERRJRSATiON LINE PANEL 


(a) e = 0.04872 ; y = 0° ; $ expanded up to e 2 ; t/C = 9.75% ; 

50 panels with cosine spacing. 

Fig. 4.14 Comparison of pressure distributions from the perturbation 
analysis method with the panel method for symmetrical 
Joukowski airfoil. 









(a) e = 0.06496 ; y = 0° ; <p expanded up to e 2 ; t/C = 13% *, 

50 panels with cosine spacing. 

Fig. 4.15 Comparison of pressure distributions from the perturbation 
analysis method with the panel method for symmetrical 
Joukowski airfoil. 
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(c) e = 0.06496 ; T = 12° ; <t> expanded up to e 2 ; t/C = 13% ; 

50 panels with cosipe spacing. 


Fig. 4.15 Concluded. 



Table 4.1 Comparison of lift coefficients of the perturbation, panel and exact 
solutions for the symmetrical Joukowskl airfoils 


Slenderness 

ratio 

e 

Exact 

C t,ex 

Perturbation 

C i»pt 

Panel 

C *,pn 

,~'_ C A»ex x 10Q% 
Jt,ex 

?***■ x ioo% 
L t,ex 




y = 6° 



0.06496 

0.71588 

0.72245 

0.67684 

0.92 

-5.45 

0.04872 

0.70234 

0.70603 

0.66472 

0.53 

-5.36 

0.03248 

0.68797 

0.68961 

0.65377 

0.24 

-4.97 

0.01624 

0.67278 

0.67319 

0.63080 

0.06 

-6.24 




Y = 12° 



0.06496 

1.42391 

1.43698 

1.34485 

0.92 

-5.55 

0.04872 

1.39698 

1.40432 

1.32128 

0.53 

-5.42 

0.03248 

1.36840 

0.37167 

1.29823 

0.24 

-5.13 

0.01624 

1.33819 

1.33901 

1.23750 

0.06 

-7.52 


o 


102 


Yt XI = .05G»iaXU/- SCRTfSs'Xm 'iPPiNr sc 
a x) =1 .ooc»xi i-x) st xi =xt i-x: *>3 on= 5 cq.o 

T/C= 3.25% fiNGLE OF flTTRCKr O.GO DEGREE 
flLPHfl= .S2'46E-03 . -.3137E-C4 3Eifi=l .OCCO 
CROSS PERnRBflriON LINE PPNEL 


-A — 
I 



(a) e = 0.02688 ; y = 0° ; <f> expanded up to e 2 ; t/C = 3.25% ; 

50 panels with cosine spacing. 

Fig. 4.16 Comparison of pressure distributions from the perturbation 
analysis method with the panel method for cambered 
Joukowski airfoil. 
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(b) e = 0.02688 ; y = 6° ; <j> expanded up to e 2 ; t/C = 3.25% ; 

50 panels with cosine sDacinq. 


Fig. 4.16 Continued. 
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(c) e = 0.02688 ; y = 12° ; <p expanded ud to e 2 ; t/C = 3.25% ; 

50 panels with cosine spacing. 


Fig. 4.16 Concluded. 
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(a) e = 0.05376 ; y = 0 ; <|> expanded up to e z ; t/C = 6.5% ; 

50 panels with cosine spacing. 

Fig. 4.17 Comparison of pressure distributions from the perturbation 
analysis method with the panel method for cambered 
Joukowski airfoil. 
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(b) e = 0.05376 ; Y = 6° ; <p expanded up to e 2 ; t/C = 6.5% ; 

50 panels with cosine spacing. 


Fig. 4.17 Continued. 


107 



(c) e = 0.05376' ; Y = 12° ; <j> expanded up to e 2 ; t/C = 6.5% ; 

50 panels with cosine spacing. 

Fig. 4.17 Concluded. 





50 panels with cosine spacing. 

Fig. 4.18 Comparison of pressure distributions from the perturbation 
analysis method with the panel method for cambered 
Joukowski airfoil. 
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(b) e = 0.08064 ; Y = 6° ; <P expanded up to £ 2 ; t/C = 9.75% ; 

50 panels with cosine spacing. 


Fig. 4.18 Continued. 




Ill 



(a) e = 0.10752' ; Y = 0° ; <j> expanded up to e 2 ; t/C = 13% ; 

50 panels with cosine spacing. 

Fig. 4.19 Comparison of pressure distributions from the perturbation 
analysis method with the panel method for cambered 
Joukowski airfoil. 









113 



X 


(c) e = 0.10752 ; Y = 12° ; <p expanded up to e 2 ; t/C = 13% ; 

50 panels with cosine spacing. 

Fig. 4.19 Concluded. 



Table 4.2 Comparison of lift coefficients of the perturbation, panel and exact 
solutions for the cambered Joukowski airfoils 


Slenderness 

ratio 

e 

Exact 

c 

i, ex 

Perturbation 

C 4,pt 

Panel 

r 

C C 

4,pt 4, ex inn* 

C - C 

4,pn 4, ex 

4,pn 

4 ,ex 

C 4,ex 

0.10752 

0.61721 

0.69115 

T = o° 
0.63346 

11.98 

2.63 

0.08064 

0.46641 

0.50658 

0.46610 

8.61 

-0.07 

0.05376 

0.31268 

0.32987 

0.30390 

5.50 

-2.81 

0.02688 

0.15689 

0.16101 

0.14812 

2.62 

-5.59 

0.10752 

1.32350 

1.40981 

Y = 6° 
1.29947 

6.52 

-1.82 

0.08064 

1.16264 

1.20984 

1.12378 

4.06 

-3.34 

0.05376 

1.99734 

1.01767 

0.95359 

2.04 

-4.39 

0.02688 

1.82841 

0.83332 

0.77878 

0.59 

-5.99 

0.10752 

2.01528 

2.11303 

Y = 12 ° 
1.94816 

4.85 

-3.33 

0.08064 

1.84614 

1.89984 

1.76670 

2.91 

-4.30 

0.05376 

1.67107 

1.69432 

1.58984 

1.39 

-4.86 

0.02688 

1.49084 

1.49649 

1.38619 

0.38 

-7.02 


100 % 
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Chapter 5 

CONCLUDING REMARKS 

The perturbation analysis method has been used to analyze the 
potential flow, due to a uniform stream, past an axisymmetric body or a 
two-dimensional airfoil. The theoretical formulations have been 
described in Chap. 2. Results obtained from the above method for 
certain classes of either axisymmetric slender bodies or two-dimensional 
thin airfoils have been compared with those generated by the panel 
method, as well as with exact solutions whenever available. 

Numerical experiments for ellipsoidal bodies and elliptic airfoils 
have shown that the higher order expansion of + In e gives better 
results, as compared with exact solutions than lower order repre- 
sentations. The perturbation analysis method for axisymmetric slender 
bodies give pressure distributions which compare well with the panel and 
exact solutions for small slenderness ratio up to e = 0.1, corresponding 
to a L/d ratio of 5. The perturbation analysis method is quite 
inexpensive, because little computer memory is required; but it is only 
applicable to slender bodies. It is suggested to be an applicable 
method as a preliminary step in modelling some simple fuselages as part 
of the analysis of a wing mounted on’ a slender axisymmetric fuselage. 

The utility of the perturbation analysis method for the two- 
dimensional airfoils is also confined to the thin airfoils with small 
slenderness ratios up to e = 0.03248, corresponding to t/C = 6.5%. The 
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perturbation analysis method seems to relieve the numerical instability 
of the panel method at some regions near the leading edge of the airfoil 
surface. The perturbation analysis method can also provide more 
accurate results for thin airfoils (i.e., e tends to a small value) as 
compared with the panel method. Moreover, the computing time and 
computer memory required are relatively small as compared with an 
existing panel method. The perturbation analysis method is recommended 
as a method which will provide reliable solutions for thin airfoils at 
low cost. Using the thin airfoil analysis as an illustrative example, 
the perturbation analysis method can be extended to analyze the three- 
dimensional wing as the next task. If the extension to the three- 
dimensional wing can be completed, the next step will be the development 
of a program for the wing design problem with the axisymmetric fuselage 
effects on a pair of wings. This study has been only a preliminary step 
for the wing -fusel age design problem. 

In addition, the perturbation analysis method can be applied to 
solve the electrostatic and magnetostatic potentials around a slender 
conducting body in a similar manner. 
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APPENDIX A 

DETERMINATION OF THE COEFFICIENTS OF a(e) and 3(e) FOR 
AX I SYMMETRIC BODIES 


Appropriate constants o k need to be determined such that 
regular at x=0. One can proceed in the following way. The 
functions g^x) are defined by 

9(x»e) = l g k (x)e 2k * {[x - o(e)] 2 + e 2 S(x)} (Al) 

k=0 K 

00 2n 

where o(s) * £ a e . Since S(0) = 0 is assumed, as x is set equal 
n=*l 

to zero in Eq. (Al), one obtains 



g k (°) 

To determine 
C* - a(e )] 2 + 

Here, a Q = - x 
a(e) 


3 a k 

explicitly. 


e S(x) = l 
k=0 


is defined. 


= I 

n=0 


2n 


, k > 1 (A2) 

one squares both sides of Eq. (Al) and has 

[j 0 9 J <X) g k-j (x >] ^ (A3) 

such that a(e) can be written as 


(A4) 


Using Eq. (A4) in Eq. (A3), the left side of Eq. (A3) becomes 
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:l [,1 1 ^] 


2k , 2^ I \ 

e + e S(x) 


Equating the coefficients of the same power of e on both sides of the 
above equation, one obtains 


k k 

l a. <v, + \ , S(x) = l g i (x)g k _.(x) 

jsQ J J j=0 J J 


(A5) 


where Is the Kronecker delta, i.e., - 1 if m=n, 5 =0 if m * n. 

mn mn mn 


For k - 0 


9 0 (x) = x 


(A6) 


For k = 1 g x (x) - -c^ + 


(A7) 


Using Eq. (A2) and with the Taylor series expansion of S(x), Eq. (A7) 
yields ai =S'(0)/4. 

For k > 2, Eq. (A5) becomes 


1 k - 1 

9k<*> * -Vs-jlj C °jVj - 9j (x)g k-j 


(x)3 


(A8a) 


Using Eq, (A2) In Eq. (A8a), one has 


1 k ’ 1 

“k ' -7 9k-j<°> 


(A8b) 


One can determine <* k recursively by using Eqs. (A8a) and (A8b) and has 


the result as follows: 
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a(£) ,M^ S '«0) e 4 




CS'(0)] 2 S ,,, (0) A S' (0)[S' • (0)3 2 
5 


e 6 + 0( e 8 ) 


(A9) 


Similiarly, g k (x) can be made regular at x=l and the coefficients 8 k are 
determined by following the same procedures as above. The resulting 
expressions are as follows: 


® 4L n O 1/2 

g(x,e) = l g k (x)e Zk * {[x - 0(e)] Z + e Z S(x)} (A10) 

k=0 * 


where S(e) * 1 - J $ n e 2n 8 0 “ 0 and 8 1 * “ S' (l)/4 
n=l 




(All) 


I Bi e + \ I s(x) = i g i (x)g k . i (x) 
j=o 3 k-j K,A j=o 3 * 3 


k > 0 


(A12) 


g 0 (x) * 1 - x 


(A13) 


?1< X > * -h * sffexT 


(A14) 


For k > 2 


3k< x > * +7pzcr jij 


(A15) 


6 k 8 j 9 k-j (1) 


(A16) 



1ZZ 


6U ) . i + e 2 . SlLUS-d), £ 4 

(A17) 

♦ ^ fcs^in^a) t s'(i)[s‘' ( i)] 2 ] e 6 + 0(£ 8, 

By using the above resulting equations, o(e) and 3(e) can be evaluated. 
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APPENDIX B 

BOUNDARY CONDITIONS IN COMPLEX PLANE 


The complex potential <fr(z) is defined in Eq. (2.35) as 
4 >(z) = «j»(x,y) + ii»(x,y) which is analytic in the z-plane outside G. The 
derivative of the complex potential with respect to z is denoted by W(z) 
and Is known as the complex velocity, i.e., 

*< z > 

_ iii 

3y “ ay 

= u(x,y) - iv(x,y) (Bl) 

where u and v are the velocity components in x and y direction 
respectively. 

From Fig. 2.2, on the surface of the airfoil, one finds that 


dy _ v 
<Jx “ u 


The equation of the profile G is given by 


(B2) 


y = e [C(x) t /S(x)] 


on 0 < x < 1 
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C'(x) 




2/S(x) 


(B3) 


Equating Eqs. (B2) and (B3), yields 


v - e u 



(B4) 


Using Eq. (Bl) in Eq. (2.39), and considering the real part of the 
resulting expression, one has the same expression as Eq. (B4). 

To satisfy the boundary condition at infinity in Eq. (2.39), one 

needs 


8 8 

Re / f(i-,e)d? * 0 and Im / f(5,e)d£ = r 
a a 

where r is the total circulation about the body. Physically, the first 
condition says that there is no fluid flow through the surface of the 
body. 
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APPENDIX C 

EXPRESSIONS FOR TWO-DIMENSIONAL AIRFOILS 

1. Linear Operators of the Integral Equation 
Linear operators LP,p=0, 1 in Eq. (2.50) are defined in Sec. 4 of 
[11], and some of the related operators are given by 

b n (x) ’ Vo ‘ a j< x > b n-j< x > ' C1 > 

- „ 22 
g(x.«) - I g„(x)c n * {[a - x - 1eC(x)r + e S(x)} (C2) 

n=0 

h(x,«) - I h (x)e n =■ {[6 - x - iC(x)] 2 + e 2 S(x)} 1/2 (C3) 

n-0 " 

h ~n,j< x > F > - fr [ (T-x-UC)-< 2k+1 ) [ F(T) (C4) 

- 2 I F(P)(X p 1 1 eCJ (T-x-UC) P j| de 
P=0 K! j]Je=0 

where T = [a(e)-8(e)] cos^0 + 3(e) and C = C(x). 

In particular, with n ■ j * 0 In Eq. (C4), yields 

ir 

h„ „(x,F) - 2 ! (Sln 2 0 - X ) -1 CF(s1n Z e) - F(x)] de (C5) 

u,u 0 

2. Determination of the Coefficients of o(e) and 3(e) 

•The functions g(x,e) and h(x,e) are defined in Eqs. (C2) and (C3) 
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where a(e) = J a e m and 0(e) - 1 - l a e m . 

m=2 m m=2 m 

Following the similar procedures mentioned in Appendix A, the general 
expressions for the determination of a(e) and 0(e) can be found as 
below: 


“o ■ 0 * «1 

«L<°> = ^ * * > 2 

g 0 (x) = x 
g x (x) = iC(x) 
g 2 (x) - - ^ 


For m > 3 


m-2 


•i lil t. 

\ ’ ' 1C ‘ (0) Vl ' 7 L, “j9‘m-j (0) 

J ^ 


(C6) 


m-2 


%<*> ■ -v ^ Cvi + Vi(«” *sr («) 


Applying Eqs. (C6) and (C7) with m=3 and m*4, one obtains 


g 3 00 - - « 3 - i 



u. / 


94 


(x) SIxL . CC(x)fs(x) . [sm f 

IfT lx 6 8x J 


3 o = ° = Bi 


h Q (x) = 1 - x 


h. (x) = - i C(x) 


\(D - 8 m , m > 2 


h 2 (x) = -0 2 + 




For m > 3 


, m-2 

e m = ‘ i8 m-l '.'(D + lj 2 S J "i-j' 1 ) 


(C8) 


V*J ■ - ♦ *0$ ‘Vi * Vi< x » 


m-2 


<■ mt- 

+ 2(l-x)' CS J ®"i-j ‘ h j (x)h m-j (x,:l 


(C9) 


Using Eqs. (C8) and (C9) with m=3 and m=4, this yields 


h 3 ( x) = - 


C(x)S(x) 
3 ZO-x) 3 


h 4 (x) = 


S 4 + S 2 


S(x) 

2(l-x) : 


[C(x)] S(x) 
2(l-x) 3 


[s(x)i; 

8(l-x) : 


By using the resulting expressions, a(e) and 6(e) can be determined. 
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APPENDIX D 

ELLIPTICAL COORDINATE SYSTEM 


Consider the equations 

x * A cosh? cosn 


y » A sinh? sinn cosg (Dl) 

z = A sinh? sinn sing 


With appropriate combination of x,y,z in Eq. (Dl), one has two families 
of elliptical (? is constant) and hyperbolic (n is constant) curves as 
shown in figure D. The fineness ratio is defined by 


thr _ length of major axis 
co * length of minor axis 


(D2) 


The system of the reference unit vectors at point P are denoted 
by e_, e , and e a , corresponding to the coordinate ?, n» $• Making 

>» T] p 

reference to [15], the scale factors can be determined 

1/2 

h * h = A(cosh 2 ? - cos 2 n) (D3a) 


h a » A sinh? sinn 

P 


(D3b) 
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Fig. D Two families of elliptical and hvoerbolic curves 
with different values of £ and n • 


13U 


The unit vectors of the desired coordinate system in terms of the 
Cartesian coordinate system are given by 

e^ = jp (T sinh? cosn + j cosh? sinn cose + Ic cosh? sinn sine) 

e = - (f cosh? sinn - 3 sinh? cosn cose - Ic sinh? cosn sine) (D4) 
n n 

e = - j sine + Ic cose 

P 


Making reference to Chap. 2 of [15], one obtains 


v . 1* - li 


(D5a) 


VXV = 


1 

h_h h. 
5 n e 


h e 5 e 

h e 
n n 

h e 5 

3 

i_ 

3 

35 

3n 

se 

\\ 

h V 
n n 

0 


(D5b) 


Using Eqs. (D3a) and (D3b) in Eqs. (D5a) and (D5b), yields 


$ + £ - «*** £ - cot, 0 

a? an 


(06) 


To find out the boundary conditions, one uses the requirements of no 
penetration and infinity conditions on Eq. (D5a), and has 


1>U»n) = constant = C' on ? = ? Q 


(07) 
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u n p 

i|>U,n) -S? — slnh 5 sin n + C" as 5 ♦ infinity (D8) 

By using Eqs. (D7) and (D8), Eq. (D6) can be solved by separation of 
variables method. 
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APPENDIX E 

SIMPLIFIED EXPRESSIONS FOR AXISYMMETRIC BODIES AND 
TWO-DIMENSIONAL AIRFOILS 

1. AXISYMMETRIC BODIES 

Using the Eqs. (3.3a)-(3.3f) with the specified S(x), one can find 
out some leading terms of f nm (x), <*(e) and 8(e) as follows: 

1.1 Ellipsoidal Body 

S(x) = 4x(l - x) 

f 10 (x) = 4 * u J l - *> 

f 2Q (x ) = 8* UJl “ &) 
f 21 (x) = - 8ir UJl - 2x) 

f 30< x) * 0 

f 31 (x) > 32, U_(l - 2 k) 
f 32 (x) - 16, U„(l - 2x) 
o(e) * z" + e 4 + 2e^ + O(e^) 
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3(e) = 1 - e 2 - e 4 - 2e® + 0(e 3 ) 

1.2 Dumbbell Shaped Body 
S(x) = 4bx(l - x)[l - bx(l - x)] , b > 2 

f 10 (x) = 4bTr UJ2x - l)[2bx(l - x) - 1] 

f 20 (x) = 4b 2 * UJC2 log(b 2 x 2 -b 2 x+b)][b 2 x(36x 4 -90x 3 +76x 2 -24x+2) 

+ b(28x 3 -42x 2 +16x-l)+2x-l] + b 2 x(168x 4 -420x 3 +356x 2 -U4x+10) 
+ b(92x 3 -138x 2 +52x-3)} 

f 21 (x) » 8b 2 ir Ujb 2 x(36x 4 - 90x 3 + 76x 2 -24x+2) 

+ b ( 28x 3 -42x 2 +16x-1 ) +2x-l ] 
o(e) = be 2 + b 2 (l + b)e 4 + 0(e 6 ) 

3(e) » 1 - be 2 - b 2 ( 1 + b)e 4 + O(e^) 

2. TWO-DIMENSIONAL AIRFOILS 
2.1 Elliptic Airfoil 

If the airfoil is symmetric and is set at zero angle of attack. 
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then y = 0, r n * 0 and C(x) = 0. Notice that the Imaginary parts of the 
functions f n (x) all vanish, which simplifies the derivation. With all 
these assumptions, one can find out f(x,e) up to 0(e ) as shown below: 

f 0 (x) = 0 

f x (x) = - 2(1 - 2x) 
f 2 (x) = - 4(1 - 2x) 
f 3 (x) = - 8(1 - 2x) 
o(e) * + e^ + 2e® + O(e^) 

0(e) = 1 - - e^ - 2e® + O(e^) 

2.2 General Joukowski Airfoil 
With the specified S(x) and C(x), all the operators can be 
evaluated, through Eqs. (3. 8a) -(3. 10c) can be evaluated and the 
constants Im f n (0) can be computed from Eq. (2.71). With this 
information, the following three f n (x) are obtained as follows: 

f Q (x) * i (2 siny)(x - 1) 

Re f.(x) = (cosy + Asiny)(x - 1)(1 - 4x) 

Im f x (x) = 2(x - l)[2(Acosy - siny)x + siny] 
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Re f 2 (x) a 0.5(x - 1){3 (cosy + Asiny)(8x 2 + 8x + 1) 

+ 2A [6(Acosy - siny)x(l - 2x) + siny(l - 4x)]} 

Im f 2 (x) = - i 0.5(x - l)[A(cosy + Asiny)(24x 2 - 16x + 1) 

+ (Acosy - siny)(24x 2 - 20x + 1) + 4siny(2x - 1)] 

Notice that f n (x) are complex -valued functions, so each f n (x) is 
represented by 


f n 00 » Re f n (x) + i Im f n (x) , n * 0, 1, 2 


Using the above results in Eq. (2.73), with the definition of lift 

coefficient C. one obtains 
A»pt 

2 3 

C 4 p£ = 2 t [siny + 0.5(Acosy + siny)e + 0.25A(cosy)e ] + 0(e ) 


The leading terms for a(e) is given In Eq. (2.67). Since this airfoil 
has a sharp trailing edge, 0(e) should be set to one. Thus, by using 
the resulting expressions, the asymptotic expansion for <|>(z;e) is given 
by 


*(z;e) = e‘ iy z 


1 

2ir 


2 

I 

k=0 


k A 
e ! 
o(e) 


ln(z-?) f k (?) d? 

C(i - c)(e - «)] l/i! 


3 

+ 0(e ) 


In particular, if A * 0 or C(x) = 0, the real part of the results are 



136 


suitable for the symmetrical Joukowski airfoil. For this case, the 
first two terms of <|>(z;e) above agree exactly with Eq. (12.9) of [10], 

3 

when his results are expanded up to 0(e ). 

The third order lift coefficient for the general Joukowski airfoil 

is formulated to compare with the results obtained from the perturbation 

analysis and the panel methods. Making reference to [15, 22, 23], the 

exact lift coefficient for the Joukowski airfoil is given by 

C # av s 8* a sirt( Y +e)/t, where y represents the angle of attack 

and a is the chord of the airfoil. The radius a, the angles 8 and 6 are 

connected by the Eq. (2.76). By applying the similar technique in Sec. 

2.1.3, a perturbation parameter Is chosen, then 

a = C(l-e 1 cos6)/(l-0.5 gjSin 2 5) and i = ^(l-ejCOSS) 2 /( l^e^cosfi) . 

By using the resulting expressions and retaining terms up to order 
2 

of , the lift coefficient can be expressed as 

C 4 ex = 2usin(Y+8)(l - E|COS$ - e 1 2 cos 2 5 - 0.5 e 1 2 sin 2 5) + O^ 3 ) 
where 8 = tan" 1 (e 1 sin&/(l-e 1 cos6)) and the lift coefficient for the 

3 

Joukowski airfoil can be computed up to 0^ ) if the angle of 
attack y» the parameters 6 and e 1 are specified. 






